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ABSTRACT 

We test third-order standard perturbation theory (SPT) as an approximation to non- 
linear cosmological structure formation. A novel approach is used to numerically cal- 
culate the three-dimensional dark matter density field using SPT from the initial con- 
ditions of two high-resolution cosmological simulations. The calculated density field is 
compared to the non-linear dark matter field of the simulations both point-by-point 
and statistically. For smoothing scales above 8 Mpc/h it shows a good agreement up 
to redshift 0. We present a simple fitting formula to relate the linear and non-linear 
density contrast that accurately recovers the non-linear time evolution for < z < 10 
at the per cent level. To address the problem of biasing between the matter field and 
the haloes identified in the simulation, we employ the Eulerian local bias model (ELB), 
including non-linear bias up to the third order. The bias parameters are obtained by 
fitting a scatter plot of halo and matter density (both from the simulation and from 
SPT). Using these bias parameters, we can reconstruct the halo density field. We find 
that this reconstruction is not able to capture all the details of the halo distribution. 
We investigate how well the large scale bias can be described by a constant and if it 
corresponds to the linear bias parameter b\ of the local bias model. We also discuss 
how well the halo-halo power spectrum and the halo-mass cross spectrum from the 
reconstructed halo density field agree with the corresponding statistics from the simu- 
lation. The results show that while SPT is an excellent approximation for the matter 
field for suitably large smoothing scales even at redshift 0, the ELB model can only 
account for some of the properties of the halo density field. 

Key words: cosmology: theory, large-scale structure, dark matter - galaxies: haloes 
- methods: N-body simulations 



1 INTRODUCTION 



fluid property can be written as the superposition of a term 
that grows with time and a second one that decays. 



Redshift surveys have revealed the existence of large-scale 
structures in the Universe: the galaxy distribution is or- 
ganized in a complex network of filaments surrounding 
underdense regions and crossing at density peaks which 
host galaxy clusters. These structures are believed to 
form through gravitational instability starting from (practi- 
cally) Gaussian fluctuations characterized by a nearly scale- 
invariant power spectrum. Mathematically, we can follow the 
growth of dark-matter density (and velocity) perturbations 
with respect to a smooth background in terms of a set of 
Eulerian fluid equations coupled with the Poisson equation. 
As long as the density and the velocity deviate only slightly 
from their unperturbed values, the evolutionary equations 
can be linearised and solved analytically. In this case, each 



At later times, however, non-linear terms in the fluid 
equations become important and it is no longer possible to 
derive an exact solution for realistic initial conditions. A 
widespread technique to compute approximate solutions is 
to use a perturbative approach and write the full solution 
as a series expansion in powers of the linear perturbation 
amplitude. This goes under the name of Stand ard Perturba- 
tion Theory (SPT, see lBernardeau et al.ll20o3 for a review). 
In linear theory, each Fourier mode of the fluid properties 
evolves independently of the others. Non-linearities in the 
dynamics correspond to couplings between modes of differ- 
ent wavelengths. When the statistical properties of the ini- 
tial fluctuations are known, a diagrammatic technique anal- 
ogous to the Feynman diagrams can be developed to com- 
pute ensemble averaged statistics. For a given statistic (e.g. 
the power spectrum), the lowest-order, non- vanishing terms 
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in the perturbative expansion give the leading expression 
while higher-order contributions provide additional correc- 
tions. Basically, this gives an expansion in powers of the 
variance of the density contrast. 

Since it does not account for multi-streaming, SPT 
should break down in the highly non-linear configurations 
that lead to the formation of collapsed structures. In hier- 
archical scenarios for structure formation this happens at 
increasingly larger scales with time. An interesting question 
is whether, and to what degree, this breakdown affects our 
ability to predict the evolution of structures on the largest 
scales. It is generally expected that SPT is meaningful in 
the so called mildly non-linear regime, i.e. for Fourier modes 
with wavenumber k < a few x 0.1 /i/Mpc. 

Comparison with N-body simulations shows that, in 
our currently favoured ACDM model and for statistics 
like the power spectrum and the bispectrum, SPT is 
rather accurate at redshifts z > 1 (for k < 0.2ft/Mpc) 
while it becomes increasingly imprecise as z — > due 
to the fact that th e variance of the density contrast ap- 
proaches unity (e.g. ICarlson. White fc Padmanabhan 2009; 
iNishimichi et all 120091 )" When this happens, all terms in 
the perturbative series are of the same importance and 
it does not make sense to truncate the expansion at 
finite order. A number of techniques have been pro- 
posed to res um or truncate the series in a more mean- 
ingful way dCrocce fc Scoccimarrol 120061 . IPietronil l200St 
iTaruva fc HiramatsrfeOOSl ). 

Another complication arises if we want to model the 
large-scale structure seen in galaxy redshift surveys, namely 
we have to account for the fact that galaxies (and dark- 
matter haloes in general) are biased tracers of the underly- 
ing mass distribution. This is clearly seen observationally as 
different galaxy types show different clustering amplitudes. 
Consistently, N-body simulations show that the distribu- 
tion of dark-matter haloes depends on their mass and also 
on other characteristics. The most common way to account 
for galaxy biasing is to assume that the relation between 
the fluctuation in galaxy counts within a characteristic vol- 
ume centred at a given location can be written as a power 
series of the corresponding volu me- averaged mass-density 
contrast l|Frv fc Gaztanagalll993h . This Eulerian local bias 
(ELB) scheme seems to give an accurate description of the 
large-scale clustering of dar k-matter haloes in simulations 
ijManera fc Gaztanagall2009l '). However, it is difficult to ac- 
curately predict the bias co efficients based on models for 
the c ollapse of fluctuations (jManera. Sheth fc Scoccimarrol 
l2010h . 

iHeavens. Matarrese fc Verd"el { 1998) have shown how 
the ELB model can be combined with SPT to compute 
the power spectrum of biased tracers of the large-scale 
structure up to next-to-leading order. Similarly, the 
dependence of higher order statistics on cosmological 
parameters and bias coefficients can be evaluated fol- 
lowin g the s ame approach (e.g. Gaztanaea fc Friemanl 
19941; iMatarrese. Verde fc Heavens! 1 19971: 



Buchalter fc Kamionkowskil 1 19991 ; ISefusatti fc Scoccimarrol 
20051 ) . 



scale used for this spatial averaging procedure the model 
does not make unique predictions. Also, the resulting spec- 
tra and multi-spectra show unwanted features at wavenum- 
bers corresponding to the smoothing scale. Renormalization 
of the bia s parameters has been proposed to a lleviate the 
problems (|Heavens et al.lll998l . lMcDonaldll2006t ). 

Irrespective of these difficulties, the ELB+SPT model 
has been used to extract the bias p arameter that best fit 
observations (e.g. IVerde et all 120021 ). Moreover, it is of- 
ten used to discuss how biasing and non-linearities mod- 
ify ba ryonic acoustic oscillations in the galaxy power spec- 
trum (jSmith. Scoccimarro fc Shethl2007l : ljeong fc Komatsul 
2006, 2007, 2009a) or to quantify scale-dependent bias- 
ing in the presence of non-Gaussia n initial conditions 
llTaruva. K ovama & Matsubara 2008; Jcong fc Komatsul 
2009bl: ISefusattil 120091; iGiannantonio fcPorcianil l201(i 



Baldauf. Seliak fc Senatorell2010l ). Many forecasts about the 



constraining power of cosmological parameters from future 
observational campaigns are based on these calculations. 

Even though each ingredient of the model has been sep- 
arately tested against numerical simulations, it still is un- 
clear what level of accuracy can be achieved by combining 
SPT with the ELB scheme to describe the distribution of bi- 
ased tracers of the cosmic mass distribution. Future galaxy 
surveys aiming at determining the origin of cosmic accel- 
eration will require models of the galaxy power spectrum 
with per cent accuracy. Can the ELB+SPT computational 
scheme satisfy such a requirement? 

In this paper, we present a direct comparison of model 
predictions against the outcome of state-of-the-art high- 
resolution N-body simulations. We follow a novel approach 
where we evaluate the SPT expansion of the mass density 
and velocity fields (up to third order) starting from the same 
realisation of the linear density field that has been used to 
generate the initial conditions of the simulation^]. This al- 
lows us to make a point-by-point comparison between the 
non-linear mass and halo overdensities while past studies 
have only focussed on two- or three-point statistics (of either 
SPT or ELB). Our analysis sheds new light on the interpre- 
tation of the bias parameters in ELB and on the effect of 
the smoothing procedure intrinsic to the ELB scheme. 

This work is organized as follows: In section[2]we present 
the principles of SPT and describe how we calculate the 
matter density field up to third order. The numerical simu- 
lations used are described briefly in section[3] In section|4]we 
compare the SPT matter density field to the simulations on 
a point-by-point basis and using one- and two-point statis- 
tics. Section [5] describes the estimation of the bias parame- 
ters, and their dependence on the halo mass and smoothing 
scales. The resulting halo density field is compared to the 
simulations, and the accuracy of the halo-halo power and 
halo-mass cross spectra is investigated in detail. We con- 
clude in section [6] 



This method has to face a technical difficulty: the ELB 
scheme only makes sense when density fields are averaged 
over a finite volume while SPT applies to unsmoothed fluc- 
tuations. Since the bias parameters depend on the actual 



1 We do not considered resummed theories fo r which the inclus ion 
of biasing has only recently been considered <Elia et alj|2010h . 
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2 STANDARD PERTURBATION THEORY 

Perturbations in the matter density field can be described 
by the density contrast at position x and conformal time r 

P(x,t) - p(t) 



5(x,t) 



P(r) 



(1) 



where p(r) is the mean density of the Universe, and the 
divergence of the velocity field v(x, r) 



?(x,t)=V-v(x,t) 



(2) 



The basic idea of SPT is that the Fourier transforms, 5(k, r) 
and #(k, r), can be written as a sum of separable functions 
of r and k: 



5(k,r) = ^A>( r )M k )> 



£(k,T) = -ff(T)^D„(r)0~ n (k), 



(3) 



(4) 



where D n (r) « [Di(r)]" w ith Di(t) being the linear growth 
factor for ACDM models |Bernardeaul 1 199481 ), and Di = 1 
today. 

In general, 8 n (k) is of n-th order in the linear density 
contrast field 5i(k): 

*»<*)=[ •• f& fe(qi),...,5i(qn) 



/ (2tt) 3 "7 (2tt) 3 
x F„(qi, . . . ,q„) <5 D (k - qi q n )| . (5) 



The kernels .Fn.(qi . . . q n ) describe the mode coupling in 
Fourier space due to the dynamical non-linearities. They 
can be calculat ed from the recurs ion relations for n > 2 
(first derived bv lGoroff et"aHll986h : 



F„(qi, . . . , q n ) = 



G m (qi 



,qmj 



[(2n+l)a(ki,k 2 ) 



^ (2n + 3)(n-l) 

) +2/3(k!,k 2 ) 

XGm-m(q?n+l, . . • , qn)] , 



G„(qi, ...,q„) = ^ 



G m (qi 



[3a(ki,k 2 



<-> (2n + 3)(n- 1) 

x F„_ m (q m+ i, . . . , q„) + 2n/3(ki, k 2 ) 

xG„_ m (q m+ i, . . . ,q„)] , (6) 

where ki = qi H + q m , k 2 = q m+ i H + q„ and Fi = 

Gi = 1. The kernels «(ki,k 2 ) and /3(ki,k 2 ) are given by 

(ki+k 2 )-ki 



a(ki,k 2 ) 



(ki,k 2 ) 



(k!+k 2 ) 2 (kl-k 2 



(7) 



In this work we use SPT to calculate the density contrast 
field up to third order. In order to meaningfully truncate the 
expansion after n = 3we need to make sure that |<5(x)| <§C 1. 
This is true for very early times and/or large scales. Statis- 
tically, the SPT density field can be described by the power 
spectra 



(27T) 3 P m „(fc)fo(k + k') = ((L(k) 5n(k')>, 



(8) 



the linear power spectrum. The next-to-leading-order cor- 
rections (also called the one-loop corrections) to the power 
spectrum are of the order of Sf, because the correlations of 
third order in Si are odd moments and vanish for the case of 
Gaussian initial conditions. This means that the following 
terms contribute at the one-loop level: 



Pl-loop(fc) = P 22 (fc) + 2Pl 3 (fc). 



(9) 



The one-loop power spectra can als o be directly calculated 
from the linear power spectru m l|Makino. Sasaki fc Sutol 
1 19921 . IJain fc Bertschingerl 1 19941 ) by inserting Eq. © into 
Eq. JSJ and making extensive use of Wick's theorem: 

P 22 (fc)= 2 J ^F 2 2 [k-c l , (l ]P 11 (\] i -c l \)P 11 ( q ), (10) 

Pi 3 (fe) = 3Pn(fc) J ^Lp n (g)it> [k.q.-q], (11) 

where F^ s ^ [k, q, — q] denotes Fg [k, q, — q] symmetrized 
w.r.t. its arguments. 

In previous studies, the one-loop corrections to the 
power spectrum were calculated using Eqs. (|10[l and (|11[) . 
In this work we explicitly calculate the third-order density 
contrast field potnt-by-potnt on a 3D grid and compare it to a 
simulated non-linear density contrast field to test the valid- 
ity of SPT at low redshifts. We start from an initial density 
field with iV 3 grid points and calculate the corresponding 
higher order density fields, which can then be evolved to 
any redshift according to Eq. J3). 

While Eq. (0 is an elegant analytical description of how 
the higher order density contrast field depends on the linear 
one, there are other formulations which are more suited for 
numerical integrations, because they do not require (n — 1) 
integrations over fe-space (the integral over the Dira c- delta - 
function can be done analytically). iMakino et al.l (|l992l ) 
showed that Eq. ([SJ can be rewritten as a recursion rela- 
tion (their Eqs. 2.15a and 2.15b), expressed schematically 
like: 



8 n — K \S\, ■ ■ ■ ,<5n-l, 



\oi, , On_l , PI, , Cti-1 

9 n = J [<5i, . . . , <5 n -i, 6\, . . . , 9 n -i 



(12) 



where (■) denotes the ensemble average and Pn(fc) is called 



where K and J are single integrals over fc-space with (dif- 
ferent) mode-coupling kernels. The starting point of this 
integration is given by the linearised continuity equation: 
6i (k) = <5i(k). This method is numerically faster because it 
takes significantly less time to calculate two single integrals 
(e.g. 82 and 62 which are needed for £3) than to perform 
a double-integral for ^3 using Eq. ([5]). This also means that 
even higher-order corrections could be calculated in the same 
amount of time and one obtains the velocity divergence field 
#(k) up to the (n — l)-th order as a side product. We find 
that this method is also more stable against numerical effects 
such as the discretization of the mode-coupling kernels. 

However, even with this hierarchical method, the inte- 
gration time still scales with the square of the number of 
grid points, because if we increase the number of grid points 
by, say, a factor of 8 (doubling the resolution in each di- 
mension), for each of these points we also have 8 additional 
grid points that contribute to each integral. This means we 
are limited by the resolution of our grid, i.e. the minimum 
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Figure 1. Left panel: <5Jj(x) (blue crosses) compared to the different SPT orders (black dotted: 5?(x), green dashed: (5j(x) + <5|(x), re d 
solid: <5g pT (x)) at redshift for the small box with R = 12 Mpc/h. Right panel: Same for the large box and R = 28 Mpc//i. (Second 
order is not shown here because it is very similar to the third order for this smoothing scale.) 



separation between grid points. The highest resolution that 
can be achieved with a common workstation is (128) 3 grid 
points. 



3 N-BODY SIMULATION 

The sim ulations that are used in thi s work are described in 
detail in IPillepich, Porciani fc Hahnl d2010t). They were run 
using the TreePM code Gadget-2 (|Springelll2005l ) in a flat 
ACDM cosmology starting from a Gaussian random field. 
The cosmological parameters are fl m = 0.279, Qb = 0.0462, 
Qa = 0.721, h = 0.701, a 8 = 0.817 and n s = 0.96. The 
output contains the 3D positions and velocities of the dark 
matter particles, which are converted into a density contrast 
field on a grid using the cloud-in-cell (CIC) algorithm. Out- 
puts for 30 time steps between z = 10 and z = logarith- 
mically spaced in (1 + z)~ x are available. The linear density 
contrast field <5i (x) , which is the basis of our calculations, is 
obtained from the initial conditions of the simulation. The 
non-linear matter field from the simulation output at z — 
we denote as 5 m (x). 

The simulations consist of two cubic volumes with 
different side lengths L, which are used to explore dif- 
ferent length scales and halo masses. The small box has 
L — 150 Mpc/h and contains (1024) 3 dark matter parti- 
cles with mass M part = 2.433 x 10 8 M Q /ft. The large box 
has L = 1200 Mpc/h and (1024) 3 particles with Mpart = 
1.246 xlO 11 M e /h. There are a total of 1,051,230 (1,953,437) 
haloes at redshift identified with the Friends-of- Friends al- 
gorithm in the small (large) simulation volume. The linking 
length is 0.2 times the mean inter-particle separation and 
each halo contains at least 100 particles. This gives a total 
mass range of 2.433 x 10 10 M e /h < M h < 1.2 x 10 15 M©//i. 
The initial redshifts are z — 70 and z = 50, respectively. 

To compare the SPT results with the simulations, we 
need to smooth the density field. We adopt a Gaussian ker- 
nel, which corresponds to a simple multiplication in Fourier 



space, and denote the smoothed fields with a superscript s: 

-(|k|ii) 2 - 

5 s (k) = e 5 <S(k). (13) 

The resolution limit of (128) 3 grid points for the SPT cal- 
culation restricts the scales that can be probed for the large 
simulation box, because each grid cell has a minimal size of 
1200/128 Mpc/h « 9.4 Mpc/h. In order for the smoothing 
kernel to extend over several grid cells, we have to choose 
R m 28 Mpc/h for the large simulation volume. For the small 
box, the resolution restriction is not significant, but we find 
that SPT does not give accurate results if 5 s (x) is of order 
unity at redshift 0. For this reason, we adopt R — 12 Mpc/h 
which corresponds to an r.m.s. density contrast of 0.6. 

We will use the output of the simulations for both dark 
matter and dark matter haloes to compare with our SPT cal- 
culation and investigate different length scales, halo masses 
and redshifts. The simulation results are of course not lim- 
ited by the SPT grid resolution, but if we want to make 
point-by-point comparisons, we have to use the same grid 
and smoothing scale to the simulated mass and halo distri- 
butions. 



4 DARK MATTER 

In the following sections, all density fields are considered to 
be at redshift 0, unless specifically stated otherwise. 

4.1 SPT Density Contrast Field 

We want to approximate the non-linear matter density con- 
trast field 51(x) with the third-order SPT field 5spt( x ) = 
5i(x) + 51 ( x ) + ^3( x )- We also have to investigate if the 
third-order expansion is actually more accurate than the 
lower orders. Fig. [T] shows the different approximations to 
<5f n (x) along a line parallel to one of the coordinate axis 
in the simulation over the box-length L (half a box-length 



© 0000 RAS, MNRAS 000, 000-000 



Testing SPT & ELB against simulations 5 



1.5 




■10 1-10 1-10 1 

<5,(x) (5 1 (x) + (5 2 (x) <5 SPT (x) 




-0.5 0.5 -0.5 0.5 -0.5 0.5 
<5,(x) <5 1 (x) + 5 a (x) <5 SPT (x) 



Figure 2. Left: Scatter plot of 5^(x) vs. linear order SPT (left panel), second-order SPT (middle panel), third-order SPT (right panel). 
Black dashed lines show S^ n /S%p T = 1- The smoothing scale is R = 12 Mpc//i. Note that we only show one out of eight points to improve 
readability. Right: Same for the large box with R = 28 Mpc//i. 



for L — 1200 Mpc/h to make the lines more distinguish- 
able). For both simulation volumes, the third-order density 
contrast gives the best approximation to the simulated mat- 
ter density contrast <5^(x). The agreement is already pretty 
good for the small box but excellent for the large box, ow- 
ing to the larger smoothing scale in the latter case. While 
the linear density contrast S\ (x) traces the overall structure 
of the non-linear field, the higher order expansion is more 
accurate. 

Fig. [2] shows a scatter plot of the non-linear field in the 
simulation and the different SPT orders for both volumes. 
Starting from the left panel, we see again that the linear 
density contrast is not a very good approximation to the 
non-linear density contrast. Although this result is not re- 
ally surprising, the plot once again shows that the linear 
approximation is especially wrong when |<T(x)| 5^ 1. The 
middle panel shows the density contrast up to second order 
and the right panel shows the density contrast up to third 
order. Note that the third-order over corrects the linear den- 
sity around 8 s w 1. From left to right the scatter around the 
black line given by <5m/<5sPT = 1 is reduced. In all three 
panels the difference to <5m( x ) is largest when |<5 s (x)| 1. 
Again, the difference in scatter between the large and small 
simulation volumes is due to the different smoothing scales 
used (12 and 28 Mpc/h). 

Fig. |3 shows the probability distribution function 
(PDF) for the linear, SPT and non-linear matter density 
contrast for the small box, with a smoothing scale R — 
12 Mpc/h. This also shows that SPT provides a good ap- 
proximation for the non-linear density contrast even at red- 
shift 0. 

4-1.1 Redshift Dependence 

So far we have only compared our SPT calculation at red- 
shift 0. Now we want to extend the comparison to higher red- 
shifts. Fig. [4] shows the conditional mean of the non-linear 
density contrast given its linear counterpart at redshift 




log(l+<5) 

Figure 3. Probability distribution function (PDF) for <5| pT (x) 
(red solid), <5J 1 ( X ) (bl ue dotted) and 5\ (x) (black dashed) at red- 
shift (small box, R = 12 Mpc/h). 

(left panel) and redshift 1.6 (right panel). The lines show the 
conditional mean density corresponding to: the linear evolu- 
tion (black short-dashed), the SPT density (red solid), the 
non-linear evolution from the simulation (blue dotted) and 
an analytic expectation from the sphe rical collapse model 
l|Bernardeaulll994bl . iMo fc Whitell 19961 green long-dashed), 
all smoothed with R — 12 Mpc/h. The hatched areas show 
the l-cr scatter contours for SPT and simulation. SPT and 
non-linear matter density agree well if <5f in is not too large, 
as expected, but the scatter around the mean SPT den- 
sity is generally larger than around 5^. The spherical col- 
lapse model (which was derived for the Einstein-de Sitter 
cosmological model) always overpredicts the non-linear den- 
sity contrast, and gets close only when <5 S < 0. The linear 
approximation overpredicts the simulation at low densities, 
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<5 lln (z=0) <5 Iln (z=1.61) 

Figure 4. Left: Mean non-linear matter density contrast from the simulation (blue dotted), SPT (red solid) and an analytic prediction 
from the spherical collapse model (green long-dashed) as a function of the linear density contrast at redshift 0, for the small box with 
R = 12 Mpc/h. The hatched areas show the 1-a scatter around the mean for 5g PT (horizontal hatching) and <5^ (45°), and the black 
short-dashed line corresponds to linear evolution. The inset shows the region around 8 B = 0, where linear theory actually overpredicts 
the non-linear evolution. Right: Same at z = 1.61. 



and underpredicts it at high densities. The right panel shows 
the same for redshift 1.61. At this redshift, SPT and non- 
linear matter density are in almost perfect agreement, but 
the scatter in SPT is still slightly larger. The spherical col- 
lapse model and the linear evolution are now closer to the 
non-linear evolution, but still show the systematics seen at 
redshift 0. 

The spherical collapse model is commonly used to com- 
pute Eulerian bias coefficients for dark matter hal oes (e.g. 
iMo fc Whitelll996l , iGiannantonio fe Porciaiill l201fj| and ref- 
erences therein), but it does not describe our data well. We 
therefore provide a new and more accurate fitting formula 
for the relation between <5^ and 8\. The redshift range we 
can study is < z < 10 in 30 bins, logarithmically spaced 
in (1 + z) _1 . For each redshift bin, we find that the param- 
eterisation 



M*) = M*) + B (z) ■ + C(a) • [Sl(z)] 2 



(14) 



gives a very good fit for all values of S\ : the ratio of the fit 
residuals to the true value is less than 5 per cent for redshift 
and less than 1 per cent for z > 6. The functions A, B and 
C depend on z in the following way: 

A{ Z ) = 0.006- °-° 54 



1 + 2 

0.045 

1 + z' 

C(z) = 0.643 - 0.012,2 



B(z) = 1.000 - 



0.122 
1 + z' 



(15) 



In the limiting case z — > oo, some of the fit coefficients do 
not show the right asymptotic behaviour (A — > 0, B — > 1, 
C — > 0) as one would expect because <5f n — > 5f, so this fit 
should not be used for redshifts outside the fitted range. 
Both density fields were smoothed with R = 12 Mpc/h. 
The choice of R does not significantly influence the fitting 
parameters, but the agreement is worse when large values of 
51 and 5^ are allowed (i.e. R is small at low redshifts). 



Another relation between the linear and non-linear den- 
sity fields that is sometimes used i s the lognormal trans- 
form ation (e.g. IColes fc Jones|[l99ll . iBernardeau fc Kofmanl 
Il995h . This has proven to yield a PDF which agrees well 
with simulations. However, this does not imply that it can be 
used as a point-by-point relation betw een the linear and non- 
linear density fields, as pointed out bv lKavo. Taruva fc Sutol 
|200ll) . We follow their parameterisation, which can be 
rewritten as 



„7(*!W-|7^ 2 ) 



with 



2 2 

7 a 



ln(l + al 



(16) 



(17) 



Here, a 2 (cr^) is t ne variance of the smoothed linear (non- 
linear) density field. We calculate a 2 and leave 7 as the fit 
parameter. This fit yields values for <5 m which are systemati- 
cally lower by « 30 per cent at redshift and still 20 per cent 
at z > 6. So indeed the lognormal transformation can not be 
used to accurately predict the evolution of an initial Gaus- 
sian field into a non-linear field. Again, this does not mean 
that the statistics of a non-linear density field can not be 
described by a lognormal field, just that one can not expect 
that, for a point-by-point comparison, the linear field used 
in the transformation corresponds to the initial conditions 
of that non-linear field. 



4.2 Matter Power Spectra 

After calculating the higher order density contrast, we now 
calculate the matter power spectra up to the 1-loop order. 
We can do this in two ways: calculating the volume average 
of the SPT density contrast (Eq. [SJ or integrating a given 
linear power spectrum Pn(fc) (Eq. I10[) , where Pn(fe) is the 
same power spectrum that was used to set up the initial 
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conditions of the simulation. This also serves as a consis- 
tency check for our grid-based calculation. The results for 
the large box are shown in Fig. [5] The small box is not shown 
because it is very similar. The lines show the results from 
Eq. (|10[) and the points show the numerical result, without 
smoothing, i.e. R — 0. The errors of the power spectra are 
calculated assuming that our fields are Gaussian: 



1000 



ol = — P 2 (k), 
riB 



(18) 



where ub is the number of modes in each fe-bin. 

The linear power spectrum from the simulation agrees 
very well with the input power spectrum until the scales 
which are affected by the resolution of the grid of the mat- 
ter field. Once we apply smoothing, this small-scale effect 
will disappear. The P22-term from the integration and the 
SPT calculation fits the simulation very well on all appli- 
cable scales. Finally, the Pi3-term from the SPT calcula- 
tion is slightly lower than the power spectrum integration 
(a factor of w 1.1). The lines that are shown in Fig. [5] are 
the integration result with limits corresponding to the scales 
which are available in our simulation box, but since the box 
only contains discrete modes and the power spectrum inte- 
gration assumes a continuous frequency spectrum, choosing 
the same upper and lower bound for both methods does 
not necessarily give the same result. A consistency check for 
both methods is if the power spectra overlap in the fc-region 
covered by both boxes. We find that the grid-based method 
gives very consistent results in this region, while there is an 
offset for the power spectrum integration. By slightly mod- 
ifying the integration limits for Pis(k), the discrepancy in 
the latter can be fixed, and both methods agree with each 
other. 

Fig. [6] shows the full third-order SPT matter power 
spectrum compared to the simulation. The lines show the 
power spectra directly calculated from the linear and non- 
linear density contrast in the simulation, which have been 
"glued together" to show both volumes at the same time. 
The points show the result from our grid-based calculation. 
For the large box, linear theory, SPT and non-linear evolu- 
tion agree until k ~ 0.1/i/Mpc, but linear theory seems 
to be the better approximation for slightly larger k. For 
the small box, the SPT matter power spectrum describes 
the non-linear evolution better than the linear power spec- 
trum for scales 0.2/i/Mpc < k < 0.3ft/Mpc, while they are 
equivalent for larger scales. We stress again that these power 
spectra (or the density contrast used to compute them) are 
not smoothed at all in Fig. [5] and Fig. [5] This explains the 
discrepancy of the SPT power spectra on small scales with 
respect to the good agreement of the density contrasts shown 
in Fig. \T\ The fact that the points deviate from each other in 
the overlap region between the large and small box is due to 
Pi3(fc), as discussed above. There is excess power on these 
scales for the full SPT power spectrum because Pia{k) < 0. 




0.01 



0.1 
k [h/Mpc] 



Figure 5. Contributions to the SPT matter power spectrum for 
the big box: Pn(fc) (blue solid line and squares), P22CO (red 
dashed line and triangles) and 2 -F13 (fe) | (green dotted line and 
stars). The points are the results from our calculations, while 
the lines are obtained from (integrating) a linear power spectrum 
(according to Eq. II 1 . 
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0.01 0.1 1 
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Figure 6. The third-order SPT power spectrum (points with 
error bars) compared to the non-linear power spectrum from the 
simulation (blue solid line) and the linear power spectrum evolved 
to redshift (black dashed line) for both boxes. Green triangles 
and red squares are used to distinguish the large and small box. 



5 HALOES 

In this section we study dark matter haloes and how they 
relate to the underlying matter density field. The haloes are 
split up in 3 mass bins for each box because the bias pa- 
rameters are very sensitive to the mass of the objects (see 



also Fig.ll2[). Table [T] shows the selection criteria for the bin- 
ning. Halo positions are assigned to the grid using the CIC 
algorithm to obtain a continuous halo density contrast field 
Sh (x) . The bins were chosen such that they are separated in 
mass, but still contain a large enough number of haloes to 
reduce the effect of shot noise. 
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Table 1. Halo mass bins: I-III for the small box and IV- VI for 
the large box. 



Bin 


Mass range [10 12 M Q /h] 


# of haloes 


I 


0.025 - 0.035 


91,511 


II 


0.5 - 2 


17,030 


111 


5 - 10 


2,693 


IV 


12.5 - 30 


376,292 


V 


60 - 150 


71,581 


VI 


300 - 1200 


7,673 



5.1 Bias Estimation 

Many different bias models exist in the literature. Here, we 
use the local Eulerian bias model 



<5 h (x)=F[<S m (x)]. 



(19) 



implying that the halo density 5h(x) depends only on the 
matter density <5 m (x) at the same point in space through 
a general function F. A choice of F which is widely used 
is that of a linear relation between <5h(x) and 5 m (x) with a 
bias factor b: 



(5 h (x) = 6<5 m (x). 



(20) 



This so-called linear biasing only changes the amplitude 
but not the shape of the power spectrum, since i\(k) = 
6 2 P m (k). It has been shown that this is not in agreement 
with the power spectra derived from numerical simulations 
(e.g. ISmith et all l2003h ■ Also, for |6| > 1, it can lead to 
<5h(x) < — 1, which of course does not make physical sense. 
However, on large scales the linear approximation is valid as 
will also be shown in section 15.2.11 

The next step would be to consider a non-linear rela- 
tion between <5h(x) and <5 m (x). If we assume that 5 m (x) is a 
smooth field, we can Taylor-expand the right hand side of 
Eq. (|19[1 around <5 m (x) as and get an expression for the 
higher-order bias coefficients bj 



<5h(x) 



3=0 J ' 



&(X) 



(21) 



where 6i is the linear bias, and bo is determined by the 
requirement (<Sh(x)) = 0. To be consistent with our third- 
order SPT matter field, we first smooth 5h(x) according to 
Eq. and truncate Eq. (|2T|) after fo 3 : 

£(x) = 6iC(x) + % (C(x)) 2 + | (C(x)) 3 . (22) 

2 D 

This means we obtain the 6j by fitting a polynomial to the 
scatter plot of <$h(x) and 5^ (x). 

However, there is a subtlety in the way the local bias 
model is applied. While the input mass density contrast field 
is smoothed on a certain scale R, and also the halo density 
field in the scatter plot is smoothed on the same scale, the 
reconstructed halo density from the fits also contains contri- 
butions from smaller scales. The reason for this is the choice 
of the filter function in Eq. (|13[) . If we used a top-hat filter 
in fc-space, this effect would vanish, because then all modes 
with k > R are set to zero (however this filter oscillates in 
real-space). But with the Gaussian filter, modes with, say, 
kR — 2 are only damped by a factor of w 0.14. When ap- 
plying the ELB model, taking the square and cube of the 



smoothed field in real-space corresponds to a convolution in 
fc-space. This unsmoothing effect can be seen especially in 
the halo power spectra (section 15.211 . and limits the scales 
where we can trust the calculated halo power spectra even 
more than the validity of the SPT approximation discussed 
in section l4~2l 

Our fitting process allows for three free parameters: bi , 
62 and &3 (60 is not a free parameter and always very close 
to so it can be neglected). We also fit the point cloud 
with fewer bias parameters and use the Akaike information 
criteriorQ (AIC) to determine the best fit without overfitting 
the data with too many parameters. We define 



5h,fj(x) EE Mm(x), 

<5 h ,f 2 (x) EE &l<5m(x) H 



^ 2 CIS , nn2 

y(M*0) • 



<5h,f 3 (x) EE ML(x) + — (CM) 



+ -g-(<M x )) ■ 



(23) 



which describes the parameterisation used for fitting and 
the notation to denote the halo density contrast Shfc (x) 
obtained from the i-th order fit. Note that 5h,f<(x) is not 
written with the superscript s due to the unsmoothing ef- 
fect discussed in the previous paragraph. The errors on the 
parameters are calculated using the jackknife method with 
8 subsamples of the density contrast fields. Scatter plots of 
5h(x) vs. <5f n (x) and the fitted polynomials are shown in Fig. 
[7] for all six mass bins. The mean and 1 — a scatter for 5^ in 
bins of 5^ are indicated by points with errorbars. One can 
see that the local bias model seems to be a rather good ap- 
proximation, although there is some scatter in the relation. 
A possible source of this could be shot noise from sampling 
discrete haloes. Assuming Poissonian shot noise, the ampli- 
tude can be estimated by calculating the number of haloes 
N s = Ph V s that are in each smoothing volume V s = 6tv 2 R 3 , 
and the mean number of haloes in that mass bin, N, from 
the third column in Table [1] 



TV 



N 



(24) 



where AN S = \/Ns (Poissonian noise) was used in the sec- 
ond step. For mass bin I with R = 12 Mpc/h, < 0.8. 
The maximum error from shot noise would be A(S^ = 0.8) w 
0.03. But one can see from Fig. [7] that the scatter does not 
increase as \5^\ increases, and that it is much larger than 
the maximum shot noise error. This means that there is an 
intrinsic scatter around th e local bias model, i.e. t he bias is 
not deterministic (see also ISomervilie et al.ll200ll l. For the 
higher mass bins with fewer haloes (right panel), the effect 
of shot noise is more prominent. There have been sugges- 
tions in the literature how to deal with stochastic ity (e.g. 
iDekel fc Laha^ll 19991 : ICai. Bernstein fc Shethll201ll ). but we 
will not include this in our modelling. The AIC values can 
be converted into a relati ve probability Wj for each model 
(see e.g. the appendix of I Porciani fc Norberd [20061 ) . with 
Wi — 1. Independent of halo masses and smoothing scale 
R, the third-order bias model always provides the best fit, 
with W3 ~ 1 while wi , W2 are lower by several orders of mag- 
nitudes. This means that although the second- and third- 
order fits in Fig. [7] differ only in regions with relatively few 



Basically a \ method which penalizes extra parameters. 
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Figure 7. Scatter plot of <5!(x) vs. ^(x). Top row: mass bins I and II (R = 12 Mpc/h), middle row: bins III (R = 12 Mpc/h) and IV 
(R = 28 Mpc/h), bottom row: bins V and VI (R = 28 Mpc/h). The lines show the polynomials of different order fitted to determine the 
bias parameters bi (black long-dashed: linear, green short-dashed: second order, solid blue: third order). The black squares with errorbars 
show the mean of in bins of <5^ and the 1-cr scatter around it as visual guidance. Note that, for readability, we only show one of every 
eight points used in the fitting procedure. 
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0.4 



log(l+<5) log(l+<5) 

Figure 8. Left: Probability distribution function for <5h( x ) (dotted blue) compared to 5j 1j / 1 (x) (black dashed) and ^^(x) (red solid). 
Right: Probability distribution function for <5h( x ) (dotted blue) compared to 5h,/ 3 ( x ) (black dashed) and 5^J(x) (red solid). Both panels 
show mass bin I with R = 12 Mpc/h. 



points, the AIC criterion clearly favours a model with 3 bias 
parameters. 

In the previous sections, we showed that <5spt ~ 5 m - So 
instead of fitting to the simulation matter density contrast, 
we can also obtain the bias parameters by fitting <5 h (x) vs. 
5| PT (x). Similar to Eq. ()23|) . we define 

5^ 2 T (x) = bi [tf(x) + 53(x)] + | (fl(x)) a , 
5^ 3 T (x) = 5i [dj(x) + 5l(x) + Jl(x)] + | [(^(x)) 2 

+ 2*!(x)*S(x)] + ^(«5 1 (x)) 3 . (25) 

J 6 

The bar is used to point out that bi and 6; are generally 
different. This is not particularly surprising because while 
we showed before that <5 m « <5|pt> the density contrast we 
fit to for the bi is truncated at the order of the fit in <5i. Con- 
sidering the third-order fit, one therefore expects the largest 
deviation between the fit parameters to occur between 63 
and 63, which is indeed the case except for the high mass 
bin V (see Table [5J). This truncation is consistent with pre- 
vious works, and using the full SPT density contrast in the 
second- and third-order fits does not improve the agreement 
between bi and bi. The errors on bi are also jackknife errors 
from 8 subsamples, and again the third-order fit is always 
preferred by the AIC. Note that the bias parameters also 
depend on the smoothing scale R. Both effects can be seen 
in Table [2] for mass bins I, III and V. We do not show the 
values of the fitted bias parameters for all possible combi- 
nations of mass bins and smoothing scales, but we comment 
on their mass dependence for a specific smoothing scale in 
section 15.1.21 The fact that the allowed range for &3 does 
not include zero in most cases is also an indication that the 
third-order bias parameter is required for the fit. 

5.1.1 Testing the Local Bias Assumption 

Now we can compare the fitted halo density contrast <5h,f 4 (x) 
to the true halo density contrast from the simulations b~i (x) . 



Table 2. Bias parameters for both third-order fits for different 
mass bins and smoothing scales R [Mpc//l]. 





Mass bin 


R 


4J vs. <5 m 


5 h vs - ^SPT 


bi 


I 


8 


0.747 ± 0.006 


0.666 ±0.008 


l>2 


1 


8 


-0.341 ±0.022 


-0.554 ±0.021 


63 


1 


8 


0.171 ±0.034 


1.542 ± 0.048 


6l 


1 


12 


0.719 ±0.009 


0.711 ± 0.010 


6 2 


I 


12 


-0.300 ± 0.034 


-0.438 ± 0.017 


63 


I 


12 


0.171 ±0.125 


0.620 ±0.111 


6l 


111 


12 


1.333 ±0.035 


1.323 ± 0.034 


b 2 


111 


12 


-0.716 ±0.140 


-1.086 ± 0.093 


bs 


111 


12 


-0.288 ±0.612 


0.066 ±0.532 


bi 


V 


28 


2.039 ±0.023 


2.037 ± 0.023 


b 2 


V 


28 


0.270 ±0.124 


-0.346 ± 0.136 


b 3 


V 


28 


-8.682 ± 1.476 


-8.556 ± 1.446 



Table 3. Correlation coefficients for the halo density contrast for 
a selection of mass bins and smoothing scales R [Mpc/h], 



Mass bin 


R 


Pair 


C 


CSp 


I 


12 


K 


5 h, fl > 


0.987 


0.989 


1 


12 




5h,f 2 > 


0.990 


0.989 


1 


12 




5 h,f 3 ) 


0.990 


0.989 


1 


12 


(Si 


*&> 


0.951 


0.961 


1 


12 


K 


C> 


0.985 


0.986 


1 


12 


(s s h 


<0 


0.989 


0.989 


111 


12 


K 


5h.fi > 


0.934 


0.956 


111 


12 




5h,f 2 > 


0.943 


0.956 


111 


12 


(Si 


<Sh,f 3 > 


0.043 


0.956 


111 


12 


K 


«&> 


0.902 


0.922 


111 


12 


(Si 


5^ > 

h,t2 ' 


0.937 


0.953 


111 


12 


(Si 


cspt ■ 


0.944 


0.957 


V 


28 


K 


*h,fl > 


0.876 


0.875 


V 


28 


(Si 


5 h,f 2 > 


0.876 


0.875 


V 


28 


(Si 




0.877 


0.875 


V 


28 


(Si 


c> 


0.866 


0.869 


V 


28 


(Si 


*h p f 2 > 


0.876 


0.875 


V 


28 


(Si 





0.877 


0.875 
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Figure 9. Probability distribution function for <5h(x) (dotted blue) compared to 5^/3 (x) (black dashed) and ^^(x) (red solid). Left: 
Medium mass bin III with R = 12 Mpc/h. Right: High mass bin V with R = 28 Mpc/h. 




Figure 10. Left: Amplitude of the halo density contrast from the simulation (blue crosses) compared to 8^ f (x). Red dotted: linear fit, 
green solid: second-order fit, black long-dashed: third-order fit. Right: Same for <5fJ^ T (x). Both panels show mass bin I with R = 12 Mpc/h. 
Note that this is the same region of the density field as in the left panel of Fig. [Tj 




50 100 150 100 200 300 400 500 600 

x [Mpc/h] x [Mpc/h] 

Figure 11. Amplitude of the halo density contrast from the simulation (blue crosses) compared to <5h f. (x). Red dotted: linear fit, green 
solid: second-order fit, black long-dashed: third-order fit. Left: Medium mass bin III with R = 12 Mpc/h. (Same region as in Fig. HOP . 
Right: High mass bin V with R = 28 Mpc/h. 
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Fig. [8] shows the probability distribution function of the fit- 
ted halo density contrast compared to the simulation. The 
left panel shows the results using a linear fit, i.e. only one 
bias parameter 61 , fitting to 8^ and S{ respectively. Neither 
fit gives a satisfactory halo density contrast, and the PDFs 
look quite different. The right panel shows the same but this 
time using the third-order fit with all three bias parameters. 
Here, both fits give PDFs which are almost indistinguish- 
able, but the agreement with the simulated halo density 
contrast is still not perfect. This discrepancy can be seen 
for all mass bins: Fig. [9] shows the PDFs of the third-order 
fits for mass bins III and V compared to the simulated den- 
sity contrast. This indicates that the local bias model does 
not capture all the properties of the halo density field. 

Fig. [10] shows again two slices of the simulation along 
one axis, extending over the box length L = 150 Mpc/ft. 
Overplotted are the fits of different order to 5^ (left panel) 
and to 5spt (right panel). In the former case, the different 
fits behave very similarly. In the latter case, the third-order 
fit (black long-dashed) seems to be closest to the simulation. 
It is clear that in both cases, the local bias model does not 
resemble the simulated halo density field on small scales. 
We also show the results for a medium and a high mass bin 
(III and V) in Fig. 111! but here we show only the fits to <5f n . 
Also for higher masses and larger smoothing scales, the ELB 
model is not satisfactory on a point-by-point level, although 
the PDFs shown before were similar. 

We also calculate the linear correlation coefficient c and 
the Spearman ranked correlation coefficient cs p which allows 
for a general, non-linear correlation between the simulated 
and fitted 5h( x ) (Table [3). The values are very close to 1, 
which shows as well that overall, the local bias assumption 
is not such a bad model for out data. From now on we will 
use the third-order fit because it is the model preferred by 
the AIC, the correlations also support this choice and it is 
consistent with the order of the SPT calculation. 



5.1.2 Mass Dependence 

Fig. [12] shows the mass dependence of the bias parameters 
(fitting 8^ vs. 5^) for the mass bins defined in Table [TJ and 
R = 12 Mpc//i. Different symbols distinguish the different 
parameters: red squares for 61 , green triangles for &2 and blue 
crosses for 63. The lines show a theoretical prediction for the 
bias parameters obtained as follows: The first step is to ap- 
ply the peak-background-split model to the halo mass func- 
ti on in the simulation (from the fitting formula presented 
in lPillepich et ai1l2010T ). The second step uses the spherical 
collapse mo del to relate the Lagrangian an d Eulerian bias 
parameters (|Giannantonio fc Porciani [2010). The two sets 
of bias parameters show the same trend with halo mass, but 
they are not i n a perfect agreemen t given the jackknife errors 
bars (see also lManera et al]|2010n . The error in 63 is strongly 
influenced by shot noise, because there are fewer points in 
the very high and low density regions which determine the 
shape of the polynomial. 

5.2 Halo Power Spectra 

In this section we test the accuracy of SPT and ELB for the 
halo power spectrum. In order to calculate the SPT halo 
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Figure 12. Mass dependence of the fitted bias parameters for 
R = 12 Mpc/h. Symbols: fit parameters, lines: predictions for the 
simulated halo mass function. 61 (red squares and solid line), 62 
(green triangles and dashed line) and (blue crosses and dotted 
line). 

power spectra, we have to re-arrange 5hj T ( x ) into terms of 
the same order in <5i(x), e.g. 

5f(k) = 6 1 5^(x) + |(^(x)) 2 . (26) 

The biased, third-order SPT halo power spectrum then con- 
sists of three terms: 

P t s 3 PT (fc) ee fl,, u (fc) + P h , 22 (£0 + 2P h , 13 (fc), (27) 

with Pb, mn (k) oc (<^ m) (k)<5< n) (k')) (as in Eq.©. In this way 
we can make sure that also the halo power spectrum does 
not contain terms of order higher than Sf. 

Fig. 1 131 shows the ratio of the reconstructed halo power 
spectrum to the one from the simulation (left panel: large 
box, right panel: small box), which has been corrected for 
(Poissonian) shot noise in the following way: 

i^fc) = itf»(fc) - ■==, (28) 

where JV is the number of halos in each bin and L is the box 
size. The top panel on each side corresponds to the power 
spectrum of 5h^, while the middle panel corresponds to 
(Jh^fJ- In the lowest panel we consider <5h,f 3 - Two effects can 
be noticed: In general, both linear and third-order SPT un- 
derestimates the power apart from the largest scales in both 
boxes. However, when kR ~ it, the unsmoothing effect sets 
in, eventually leading to a diverging ratio. The lowest panel 
shows that this behaviour is not caused by SPT, because 
dhjj an< i <5h,f 3 have very similar power spectra. 

We also tried to reverse the order of smoothing and 
fitting in the following way: Determining the bias parameters 
from a scatter plot of the unsmoothecjf] halo and matter 

3 That is without applying Eq. (15) . Of course, the fields are 
© 0000 RAS, MNRAS 000, 000-000 



Testing SPT & ELB against simulations 13 



PL, 



PL, 



1.5 
1 

0.5 

1.5 
1 

0.5 

1.5 
1 

0.5 



n — i — i — i I i i i i. 



"J J jjj— $ — I— ffl- 



H 1 — i — i m i n i m: 



\ if — I gj i 



"1 — m — s — w 



H 1 — i — i i i i I 1 1 1 1 : 



-J if — j — | — j — i- 



_L 



_l I I I I I I I 



> 1 

^ 0.5 

* 2 

^ 1.5 

L- i 

fc 0.5 

« 1.5 
eu 

>. 1 

57 0.5 



~ i — i — i — i — i — i — i — i — i i i i i i i i 1 1 1 1 1 1 1 1 1 1 iji 



H — I I I I I I I I I I I I I I 1 1 1 Mi l I P 



H — I I I I I I I I I I I I I I I 1 1 1 1 1 1 1 1 1 I I 



_i i i i i i iii 



1 1 1 1 1 A 



0.01 



0.05 
k [h/Mpc] 



0.1 



0.1 



0.2 
k [h/Mpc] 



0.3 



Figure 13. Ratios of halo power spectra from our fits to the simulation. Left panel: big simulation volume (mass bin IV, R = 28 Mpc/h), 
right panel: small simulation volume (mass bin I, R = 12 Mpc/ft). From top to bottom: linear fit to the linear matter density; third-order 
fit to full SPT density; third-order fit to the non-linear matter density field. 
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Figure 14. Ratios of halo-matter cross spectra from our fits to the simulation. Left panel: big simulation volume (mass bin IV, 
R = 28 Mpc/h), right panel: small simulation volume (mass bin I, R = 12 Mpc/h). From top to bottom: linear fit to the linear matter 
density; third-order fit to full SPT density; third-order fit to the non-linear matter density field. 



densities, then multiplying the (unsmoothed) matter density 
with these new bias parameters and applying Eq. (|13|l to the 
resulting biased halo field. The halo power spectra from this 
fit do not show the unsmoothing effect by construction, but 
are also incompatible with the simulation halo power spectra 
on any scale. In fact, not even the linear bias parameter from 
this new method agrees with the previous estimates, so that 



always smoothed at least on the scale of the grid used in the CIC 
algorithm. 



the halo power spectrum is different from the simulation 
even on large scales where the smoothing should have no 
effect (because Ph « b\P m ). We therefore conclude that (not 
surprisingly) the order of smoothing and fitting can not be 
reversed, and one has to live with the limitations of the non- 
linear local bias model that re-introduces some small-scale 
fluctuations. 

We show the ratio of the halo-matter cross spectra in 
Fig. 1141 the order of the panels is the same as in the previous 
figure. The cross spectrum is not affected by shot noise, and 
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Figure 15. Left: Halo power spectrum and residuals for <5h f 3 (x) (upper panel, mass bin I, R = 12 Mpc/h). Dashed blue line: P^(k), 
solid red line: Ph,f 3 [k), dotted black line: PaaCOi points: the cross-term 2 Pa s h fg (A)- Open squares show positive values, solid triangles 
show negative values of the cross-term. Ratio between Ph,i 3 (k) an d Ph(k) (lower panel). Right: Same for the halo matter cross spectrum. 
Dashed blue line: PhmCO, points: i"A<5 m (k), solid red line: Pg h f f m (fc). Note the change in sign at k as 0.3 h/Mpc. 



ELB spectra agree with the simulation over a larger range of 
scales for the big box. Contrary to the halo power spectra, 
the unsmoothing effect leads to an additional loss in power 
for the cross spectra on the very small scales, particularly 
evident for the small box. Note that the best agreement 
between simulation and fits is actually achieved using the 
linear bias parameter and the linear matter power spectrum 
(top panel). For the large box, this is not suprising because 
the scales in consideration are large enough for linear the- 
ory to be a good approximation. For the small box, linear 
theory is underestimating the power for the small scales, 
as seen before in Figs. [6] and 1131 but only on scales where 
the cross spectra from the higher-order fits are also too low. 
This means that the unsmoothing effect (which does not 
affect the top panel) has an even stronger influence on the 
halo power spectra from SPT+ELB than the deviation from 
linear theory on small scales, at least within the accuracy of 
our simulations. 

In order to investigate the systematic effects seen in 
Figs. \T3\ and 1141 we now consider the influence of the fit 
residuals on auto and cross spectra, which include both the 
intrinsic scatter as well as the unsmoothing effect from the 
local bias scheme. We define the residual A(x) for the third- 
order fit by 

A(x) = <%(x)-5 hl f 3 (x), (29) 

from which 

Ph(fe) = Ph.i 3 0) + 2 P A <f M3 (k) + PaaO), 
Phm(fc) = Ps Kt3 s m (k) + PAS m (k), (30) 

follows for the halo power spectrum and the halo matter 
cross spectrum. The left panel of Fig. [15] shows the dif- 
ferent terms for Ph(k): Ph,f 3 (k) (solid red) is very close to 
Ph(k) (dashed blue) until k « 0.3 h/Mpc. The cross-term 



2 Pa <s h f (k) (green points) can have both positive and neg- 
ative values. This is indicated in the figure by the differ- 
ent symbols: open squares for positive values, solid trian- 
gles for negative values. Finally, Paa(&) is shown with a 
black dotted line. Generally, 2 Pa s h f (k) and Paa (k) have 
a small amplitude, but their contribution becomes impor- 
tant on smaller scales. The lower panel shows the ratio of 
the fit and the simulated halo power spectrum for compari- 
son. The region where Ph,f 3 (k) is systematically lower than 
Ph(k) roughly corresponds to the region where the cross- 
term is positive. 

The right panel of Fig. [15] shows the terms contributing 
to Phm(fc) (blue dashed): the cross spectrum of the fit and 
the dark matter (red solid) and the cross spectrum of the 
residuals and the dark matter (black points, squares positive 
values, triangles negative values). In contrast to the fitted 
halo power spectrum which rises on small scales, the fitted 
cross spectrum becomes negative at around the same scales, 
indicated by the "hole" in the red solid line (for smaller 
scales we then show —Ps b f <5 m ). The lower panel shows the 
ratio of the fit and the simulation cross spectrum. As be- 
fore, squares correspond to positive values of this ratio, and 
triangles show where it is negative. 

We conclude that the halo auto and cross spectra com- 
puted with the ELB deviate from the simulation on scales 
where kR ~ n. This affects especially the auto spectra, and 
is not related to using SPT instead of the non-linear matter 
field. 

5.2.1 Biased SPT on Large Scales 

iHeavens et alj (1998) discuss two effects on the halo power 
spectrum that are caused by using SPT. The basis for their 
analysis is as follows: as for the case of the one-loop mat- 
ter power spectra, one can also express the SPT halo power 
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Table 4. Comparing &i, 61 , 6 c ff , fe^j and for different mass bins and smoothing scales R. (■) here denotes the average over a specific 
fc-range (see text). R is given in Mpc/h. Errors are jackknife errors using 8 subsamples. 



Mass bin R 61 61 (6 eff > (b^ B ) (6^) 



I 


6 


0.773 


± 


0.005 


0.573 


± 


0.011 


0.684 ± 


0.051 


0.606 


± 


0.040 


0.623 ± 


0.028 


11 


6 


0.928 


± 


0.011 


0.680 


± 


0.013 


0.756 ± 


0.057 


0.484 


± 


0.067 


0.800 ± 


0.032 


111 


6 


1.476 


± 


0.030 


1.108 


± 


0.043 


1.198 ± 


0.091 


0.722 


± 


0.112 


1.100 ± 


0.061 


I 


12 


0.719 


± 


0.009 


0.694 


± 


0.001 


0.684 ± 


0.051 


0.545 


± 


0.031 


0.703 ± 


0.044 


IV 


28 


1.279 


± 


0.008 


1.278 


± 


0.001 


1.253 ± 


0.076 


1.056 


± 


0.042 


1.243 ± 


0.059 


IV 


50 


1.244 


± 


0.028 


1.241 


± 


0.028 


1.253 ± 


0.076 


1.185 


± 


0.060 


1.237 ± 


0.054 


V 


28 


2.040 


± 


0.023 


2.037 


± 


0.022 


1.958 ± 


0.119 


1.553 


± 


0.295 


1.947 ± 


0.095 


V 


50 


2.029 


± 


0.042 


2.016 


± 


0.039 


1.958 ± 


0.119 


1.697 


± 


0.111 


1.972 ± 


0.089 


VI 


28 


3.704 


± 


0.054 


3.615 


± 


0.049 


3.520 ± 


0.222 


4.739 


± 


0.134 


3.662 ± 


0.159 


VI 


50 


3.761 


± 


0.088 


3.739 


± 


0.087 


3.520 ± 


0.222 


4.461 


± 


0.296 


3.877 ± 


0.155 



spectra as integrals over products of th e linear matter power 
spect rum and the bias parameters (|Jain fc Bertschingerl 
1 1994 ) • The Ph,22 and Ph,i3 terms defined before can then 
be written as 



Ph, 22 (fe) 



d 3 q 
(2tt)3 



Pn^Pndk-qi; 



^'(q.k-qj + l 



+| + b 2 p 2 (s) ( 



-q,k) 



(31) 



and Ph,n(fc) = b\P\i{k). By taking the limit k — > of the 
foregoing equations, one can study the behaviour of the SPT 
halo power spectrum on large scales: 

(i) Is the large-scale bias not fei? W hile the linear 
bias m odel should be valid on large scales, iHeavens et al.l 
(1998) predict that the large-scale bias is not 61, but can be 
approximated by an effective bias 



b e a = \ bf + bi 



b2 + 63 



(32) 



where a R is the variance of the smoothed linear density field 
at redshift 0: 



2 

Or 



d 3 g 



Pn(g)e 



-(gR) 2 



(33) 



Eq. (|32[) tells us that the difference between bfg and 61 de- 
pends on the choice of the smoothing scale through a\. How- 
ever, 61 is supposed to describe the large-scale behaviour of 
the power spectrum, where the smoothing scale should not 
have any effect. The value of a\ will be large for smaller R 
and fe 2 is negative in many cases, so the term in parenthesis 
can even become negative, leading to 6^ ff < 0. (Although 
this never happens for the smoothing scales we investigate.) 
We can compare this with the effective bias directly from 
the simulation power spectra 

(34) 

on large scales (k < 0.02 /i/Mpc for the large box and 
k < 0.15 fo/Mpc for the small box), where this ratio is found 
to be constant within the errorbars. Note that Ph has been 



corrected for shot noise, and for this specific ratio, the de- 
pendence on the smoothing scale cancels out. Examples of 
the different bias parameters are shown in Table [4] for mass 
bins in both boxes and several smoothing scales. Here, bi 
and bi are the linear bias parameters from the third-order 
fits, and bf s was also calculated from these values. We al- 
ways find the large-scale bias to be very close to 61 for rea- 
sonably smoothing scales. Note that the value for 6 c ff for 
the small box (mass bins I-III) is only an approximation for 
the large-scale bias, because the box is so small. Compar- 
ing 61, bi and (6 e ff) f° r the low mass bins, it is clear that 
the different estimates do not agree with each other if the 
smoothing scale gets too small (R — 6 Mpc//i). Note also 
that exchanging Ph with one of the fitted halo power spec- 
tra in Eq. ()34[) does not significantly affect the value of b e ff 
(exchanging P m — > Pspt for the SPT fit as well). 

Extending the comparison to b^ s reveals a large discrep- 
ancy to the previous values over the whole mass range and 
different smoothing scales. Even though the errors on 6^ are 
rather large (owing to the large errors in 63), the estimates 
are quite different from the fit parameters and 6 c ff. 

In concordance w ith the derivation of Eq. (|32|l in 
IHeavens et alj l|l998h . we define an additional effective bias 
for the SPT fit: 



err 



Pll ' 



(35) 



which is also averaged over the same scales as 6 e ff- The 
difference to Eq. (|34|l comes because the denominator con- 
tains only the linear power spectrum. However, the values 
of b s J^ do not differ much from b e g because on large scales 
P m « Pspt « P11. 

We conclude that b e s matches the fit parameters bi and 
61 within the errorbars. Eq. (|32p was derived neglecting the 
contribution of Ph, 22 (A:), which seems not to be a valid ap- 
proximation, as the values for b^g do not agree with 6 e ff and 
. It is also important to note that the agreement of the 
different bias estimations gets worse when smaller smooth- 
ing scales are considered. This implies that using the bias 
parameters to infer halo masses (by choosing R to corre- 
spond to the Lagrangian radius of the halo) can lead to 
wrong results. 



4 Meaning that the SPT assumptions about the smallness of 5 
are still valid at z = 0, requiring R > 8Mpc//i. 
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Figure 16. Left: Different contributions to 22W (see text) and their sum (blue solid) for mass bin IV with R = 28 Mpc//i. The 

fa 2 ~ b 

constant power on large scales is dominated by {S\ (k)<5|(k')) (green line). Right: The contribution A ■ <5^(x) (solid red, with A = -y- ■ 5 
for enhanced visibility), which coincides with the peaks and troughs of the underlying SPT matter density field (dashed blue) as suggested 
by Heavens et al. (1998). 



(ii) Is the SPT halo power spectrum constant on 
very large scales? From Eq. (13ip we see that the Ph,22(k)- 
term is not directly proportional to Pn(fc) which falls off 
as k — s> 0. It follows that this term will eventually domi- 
nate, and it can be shown to lead to a constant halo power 
spectrum on very large scales. We can study this behaviour 
using the SPT halo density contrast from the third-order 
fit. Even with our large simulation volume, we can only see 
that Ph,22(fc) tends to a constant, but the scales where it 
actually dominates are out of reach. However, we can look 
in more detail at the Ph ; 22(fc)-term, to find out why it be- 
comes constant. From Eq. ((26}, we can see that there a 
three terms that contribute, which are shown in the left 
panel of Fig. 1161 &fP22(fc) (red dot-dashed), the cross term 
bib 2 {5 s 2 (k)6f(k')) (black dotted) and ^{df(k)Sf{k')} (green 
dashed). The latter ter m clearly dominates P h,22(fc) (blue 
solid line) on large scales. Irleavens et alj (| 19981 ) suggest that 
this constant power comes from the peaks and troughs of the 
underlying density field, which we can confirm: The right 
panel of Fig. [16] shows the SPT matter density contrast 
<5§ PT (x) and 5f(x) along a line in the simulation volume. 
The latter term can have either sign depending on 62, here 
we show the mass bin IV where 62 < 0. Note that ^ 8\ has 
been multiplied by a factor of 5 to make the effect more 
visible. 



6 CONCLUSION 

In this work we have followed a novel approach by evalu- 
ating the SPT expansion of the mass density and velocity 
fields up to third order, starting from the same realisations 
of the linear density field that has been used to generate the 
initial conditions of two N-body simulations. This allowed 
us to make a point-by-point comparison between the non- 
linear mass and halo overdensities (using the Eulerian local 
bias model), while past studies have only focussed on two- 



or three-point statistics. Our results can be summarized as 
follows: 

• We found that SPT provides a good approximation to 
the density field up to redshift 0, for smoothing scales R > 
8 Mpc/h. This ensures that the linear density contrast is 
typically less than one. 

• We have compared the redshift evolution of the non- 
linear matter density field of the simulations and SPT with 
the prediction of linear theory and the spherical collapse 
model. We found that SPT is very close to the simulated 
density for all redshifts, while both linear theory and the 
spherical collapse model are in po or agreement with the sim- 
ulations. The lognormal model bv lKavo et al.l l|200ll ) gives a 
good PDF for the non-linear density contrast but fails when 
used for a point-by-point comparison. 

• In Eq. (|14|l . we presented a simple fitting formula for 
the non-linear density contrast as a function of the linear 
density contrast, which is accurate at the per cent level over 
the full range of redshifts available to us (0 < z < 10). 

• We compared the SPT matter power spectra with the 
linear and non-linear matter power spectra from the simula- 
tion. On large scales, linear theory provides a good approx- 
imation to the simulated matter power spectra, but SPT is 
superior on smaller scales, up to the maximum wavenumber 
we can probe, k « 0.3 ft/Mpc at redshift 0. 

• Assuming a deterministic, Eulerian local bias model 
(ELB) with up to 3 free parameters, we obtained values 
for these parameters by fitting polynomials to a scatter plot 
of the smoothed matter and halo density contrast. We find 
that the third-order bias model is always preferred by the 
data over models with less parameters. The reconstructed 
halo density is similar to the simulation, but the ELB can 
not accurately reproduce the simulated field. 

• We found that the bias parameters from fitting the halo 
distribution to the simulated matter density field and to 
the SPT density field generally differ, but the correspond- 
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ing power spectra are very similar. The mass dependence 
of the bias parameters shows the same trend as theoretical 
predictions based on the peak-background-split approach. 

• We compared both the halo-halo power spectrum and 
the halo-matter cross spectrum from the fits with the ones 
from the simulation. The cross spectrum is in better agree- 
ment with the simulation than the auto spectrum which 
deviates from the simulation on scales much larger than the 
smoothing radius. This is not related to using SPT instead 
of the non-linear matter field. 

• We have investigated two effects on the halo power spec- 
tra, which have been predicted by analytic considerations 
of SPT and the local bias model. First, we estimated the 
large-scale bias 6 e ff using the halo and mass power spectra 
from both the simulation and the fits. This large-scale bias 
is compatible with the linear bias parameter obtained from 
the polynomial fit if we make sure that the SPT assumptions 
are not violated, choosing the smoothing radius R such that 
5 <5cC 1. This suggests that the effective bias does not re- 
quire perturbative corrections, contrary to previous results 
based on SPT. Second, we determined the origin of the con- 
stant shot-noise term on very large scales, which is caused 

by ^-(SfSl) as predicted in lHeavens et all jl998fh 

In summary, our study shows that SPT is a suitable ap- 
proximation for the matter field even at redshift 0, provided 
a large enough smoothing radius is adopted. However, the 
Eulerian local bias model can not fully describe the halo 
density field, which is most evident from our point-by-point 
comparison in Fig. 1101 
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